function [beta_new,r] = newton_raphson_round(epsNR,maxNR,n,m,G,Gprime,Gprimeprime,...
        weight_function,sup_norm,beta_old,y,v,round_digits)
    
% Modified Newton-Raphson algorithm to round to a certain number of digits
% throughout the calculations

% specified number of digits each time
%
% Sanjeeva Balasuriya
% May 2024
%


% Function which performs the Newton-Raphson iteration to find beta
% Accepts either weighting function (good or bad)

% Initialize
%
p(1:n) = NaN;
deriv_sum = zeros(m,m);
r=1;
beta_diff = 1e10;

% Newton-Raphson iteration
%
% Perform while number of iterations is less than a specified value
% (maxNR), and the difference in the calculated betas from the previous two
% steps is greater than a specified value (epsNR)
while r<=maxNR && beta_diff>epsNR
    vector_sum = zeros(m,1);
    
    % Compute the vector term
    for i = 1:n
        p(i) = round(G(beta_old * v(i,:)'),round_digits);
        vector_sum = vector_sum + round(Gprime(p(i)) * v(i,:)' ...
            * weight_function(p(i),y(i)),round_digits);
    end
    
    % Compute the matrix (derivative) term
    for kk = 1: m
        for jj = 1: m
            deriv_sum(kk,jj) = 0;
            for i = 1: n
                wp = round(weight_function(p(i),y(i)),round_digits);
                Gp = round(Gprime(p(i)),round_digits);
                deriv_sum(kk,jj) = deriv_sum(kk,jj) + ...
                    round((wp * Gprimeprime(p(i)) - wp^2 * Gp) ...
                    * Gp * v(i,kk)*v(i,jj),round_digits);
            end
        end
    end
    
    % Invert matrix
    inverse_deriv = round(inv(deriv_sum),round_digits);
    
    % Newton-Raphson step
    beta_new = beta_old - transpose(round(inverse_deriv * vector_sum,round_digits));
    
    % Stop if already diverged
    if(sum(isnan(beta_new))>0)
        r = NaN;
    end
    
    % Prepare for next step
    beta_diff = round(sup_norm(beta_new,beta_old),round_digits);
    beta_old = beta_new;
    r = r+1;
end
